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Abstract 

ThermoAcoustic Tomography (TAT) is a promising, non invasive, 
medical imaging technique whose inverse problem can be formulated 
as an initial condition reconstruction. In this paper, we introduce a 
new algorithm originally designed to correct the state of an evolution 
model, the back and forth nudging (BFN), for the TAT inverse problem. 
We show that the flexibility of this algorithm enables to consider a 
quite general framework for TAT. The backward nudging algorithm is 
studied and a proof of the geometrical convergence rate of the BFN is 
given. A method based on Conjugate Gradient (CG) is also introduced. 
Finally, numerical experiments validate the theoretical results with a 
better BFN convergence rate for more realistic setups and a comparison 
is established between BFN, CG and a usual inversion method. 

1 Introduction 

ThermoAcoustic Tomography (TAT) is a hybrid imaging technique that uses 
ultrasound waves produced by a body submitted to a radiofrequency pulse, 
uniformly deposited throughout the body. The absorption of this initial en- 
ergy causes a non-uniform thermal expansion, leading to the propagation of 
a pressure wave outside the body to investigate. This wave is then measured 
all around the body with piezoelectric transducers or, more recently, thanks 
to interferometry techniques. 

It appears that the absorption of the initial pulse is highly related to the 
physiological properties of the tissue (l9| . As a result, the magnitude of the 
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ultrasonic emission (i.e. thermoacoustic signal), which is proportional to the 
local energ y d eposition, reveals physiologically specific absorption contrast. 



See [38j, [2^, |3j, |4j| for an introduction to the experimental setup. Considering 
that the initial illumination is a Dirac distribution in time, the problem of 
recovering the absorptivity of the investigated body from the thermoacoustic 
signal is equivalent to recovering the initial condition of a Cauchy problem 
involving the wave equation from the knowledge of the solution on a surface 
surrounding the imaging object [ST^]- Note that the experimental constraints 
do not allow the acquisition surface to completely surround the tissues to 
investigate, so that one cannot expect a better situation than measurements 
on a half-sphere (as in breast cancer detection for example). Moreover, 
several experimentation constraints limit the number of sensors, so that the 
methods of reconstruction considered must be robust to the space sampling. 

Man y w orks dealing with TAT have been achieved during the last decade 
(see [23!| for an overview, or 

and many of the methods used in these 
papers strongly depend on a set of additional assumptions, among which: 



• the homogeneity of the tissues, leading to a constant speed of sound 
for the pressure wave; 

• the lack of frequency- dependent attenuation, so that the pressure wave 
obeys the classical, undamped, wave equation; 

• the complete data situation, where the acquisition surface is considered 
to enclose the whole body. 



The study of TAT in the case of an homogeneous medium without atten- 
uation gave rise to explicit inversion formulas. These latter usually require 
that the acquisition set is a closed surface, even if they can be approximated 
in the limited view situation (where the surface is a half-sphere for example) . 
Among these, filtered backprojection formulas, even though they have been 
very successfully implemented in several situations, seem to face some issues 
when the source has support partly outside the observation surface, or when 
the data sampling is not sufficient (see (l^ . E^). Fourier's type formulas, 
which offer a much better numerical efficiency, apply on very specific ge- 

or imply some additional approximations, like an 



ometries (see 
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interpolation from spherical to cartesian coordinates P, U, [13, 2A, 45, 47 1. 
Even though some of these techniques give rise to very efficient reconstruc- 
tion schemes (especially in [24]). so far it is not clear if they can be extended 
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to less restrictive acquisition geometries or to an attenuated wave with non 



constant speed (see [25| for an answer about the acquisition surface). 

Moreover, since the work of Burgholzer et al (see 0|), the time reversal 
method (first suggested in (l^ ) has been apphed with success to the TAT (see 



e.g. [l6j, [iTjl). Even though the first known results on this method required 
a complete data framework, recent works by Stefanov et al extended it to a 
Neumann Series method in a fairly general situation (incomplete data, known 
variable sound speed and external source) with very good results in 3^, 39|. 

In this paper, we show that the framework of TAT lends itself to the appli- 
cation of some data assimilation techniques. These latter, mainly used in 
geophysics, aim at correcting the state of an evolutionary model by means 
of data, in order to obtain a good approximation of the real state (see (4o| 
for an introduction) . The particular data assimilation method proposed here 
is based on a nudging technique: given an evolution model of the state and 
direct observation (our data), it consists in adding, inside the model equa- 
tion, a newtonian recall of the state solution to the observations (or data), 
which is usually called the feedback or nudging term As we shall see, 
from a practical point of view, this method can be successfully used to man- 
age the usual issues of the TAT inverse problem as incomplete data, external 
source and variable sound speed (when given, however). So far, however, the 
theoretical convergence result for the nudging technique is based on a clas- 
sical result about stabilization of the wave equation 2^, 3^, which requires 



somehow a geometric optics condition (see [6] or [27| for an overview). This 
geometric condition, whose verification is not a simple matter, is not auto- 
matically satisfied in an incomplete data framework or with a variable sound 
speed. Consequently the proof provided in this article only stands in a very 
favorable situation (3d case, constant speed and complete data), but does 
not depend on usual stabilization techniques and yields explicit convergence 
rates. 

Different nudging terms can be found in the literature, and most of them 
have been used to assimilate data in physical oceanography as in 21, 2^, 30 1 



and more recently in Fundamental articles offer the basis of most popu- 
lar techniques as Kalman filter | |2Q] . Kalman-Bucy filter Q and Luenberger 
observer 



Since the basic nudging method showed its limits in real conditions by using 
the data once, Auroux and Blum proposed to extend the method by adding 
a resolution backward in time and iterating the process, thus defining the 
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back and forth nudging (BFN) algorithm. They presented this technique 
in 0, 0| and implemented it with computationally efficient results in com- 
parison with traditional data assimilation methods, like optimal filters (from 
extended Kalman's to particle filters) or other optimal minimization tech- 
niques (variational methods), which appear costfull respectively in mem ory 
needs and computing time - see the brief introductions of Talagrand in [40|, 



41l . I42l | for further details about data assimilation and a relevant bibliogra- 
phy. Recently, a more general formulation of the BFN using observers has 
been presented by Ramdani, Tucsnak and Weiss in js^l- 

As inverse problems are often solved by means of variational techniques, 
we also introduce a least squares method based on the Conjugate Gradient 
(CG). 

The two methods introduced are compared to the Neumann Series presented 
in |33l |. which is one of the best existing method which fits the variable speed 



case and partial data settings. 

This paper is organized as follows. In Section [21 we introduce the TAT 
inverse problem in a quite general form, i.e. without the use of the usual 
assumptions, then we describe the BFN algorithm and we state the main 
result of this article: the convergence Theorem |3l Then Section [3] is devoted 
to the proof of this result. In Section [H we describe the variational method 
that uses CG algorithm. Finally, we give numerical results in Section [5l 



Notation 

In the following, we shall use this notation: 

• The open ball with center x G E,^ and radius r > is denoted 
by B{x,r); 

• Its boundary, the sphere with center x and radius r, is noted S{x,r); 

• For every x G and < ri < r2, we define the spherical shell of 
center x and radii ri, r2: 

A{x,ri,r2) = {y G E^|ri < \\x - y\\ < , 

and if X = 0: 

A{ri,r2) = A(0,ri,r2); 
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• For any set S of and ri > 0, we denote by T{S,ri) the set S 
extended to the points at a distance lower than ri from S: 

T{S,ri) = {y £ ]R^|3x G 5, ||y-x|| < ri} . 



2 Presentation of the method 
2.1 The general TAT problem 

In the foUowing, we will denote by Pexact(a;,i) the pressure wave resulting 
from the thermal expansion of the body. We will make the assumption 
that the measurement process of this pressure is subject to some perturba- 
tion Pnoisej such that the actual data can be written Pdata ■= Pexact+Pnoise- In 
thermoacoustic tomography, the data are nothing but {pdata{x,t)\x G 
where ^5^ is a surface surrounding the body to investigate. The pressure Pexact 



satisfies the Cauchy problem (see [37| for a detailed calculation): 



Lpc^a.ct{x,t) = Edepix,t), {x,t) G X [-l,Oo), 
Pexa.ct{x, -1) = 0, X G R3, 

dtPcxnctix, -1) = 0, X G R^, 

where L is the differential operator governing the wave, most likely the wave 
operator, or a damped wave operator, and -E^ep is the energy deposited in 
the body around the time t = (so we have to consider a negative initial 
time, e.g. t = —1). This energy can be written: 

dj 

Edep = e{x) — {t), 
/3c 

e{x) := — Iem{x)'^{x), 

Cp 

where e is called the normalized energy deposition function. Here /3 is the 
thermal expansion coefficient, Cp is the specific heat capacity, c the sound 
speed and lem is the radiation intensity of the initial energy pulse. All these 
parameters, including the time shape j of the pulse, are known and we see 
that the knowledge of E^ep is sufficient to compute the absorption density ^ 
inside the body, which is our purpose. 

So far, the general thermoacoustic problem can be formulated in the following 
way: 
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Let Pexa.ctix,t) be solution of: 

dj 

i^Pcxact(x,t) = fo{x) — {t), {X,t) G ]R3 X [-l,oo), 

Pexact(a:,-1) = 0, x£B?, (1) 

dtPexi,ct{x, -1) = 0, X G E^, 

where j is known and /o is the object to reconstruct, with com- 
pact support in the unit ball -6(0, 1). Can we compute /q, or a good 
approximation of /o, from the knowledge of pdata = Pexact +Pnoise 
on a surface surrounding -6(0, 1) ? 



2.2 Reduction to the homogeneous case 

As we have seen, even though it is not possible in practice, the initial energy 
deposition is meant to be a Dirac pulse in time at i = 0, that \s ^ = 5q. 
Let be solution, which is no longer to be considered for negative time 
values, where it is the null solution, of the Cauchy problem: 

L(?o(a;,t) = /o(x)5^(t), {x,t) G E3 X R+, 
go(x,0)=0, xGR3, 

which is equivalent to: 

LqQ{x, t) = 0, (x, t) G X ]R+, 

'?o(x,0) = /o(x), XGR3, (2) 

dtqo{x,Q)=Q, xGE^, 

according to Duhamel's principle (see provided that L is a (damped) 

wave operator. Then Pdata is solution of if and only if pdata = j*tQo, which 
means that a deconvolution operation on the thermoacoustic signal Pexact 
leads to the knowledge of the solution qq of ([2]) on the same surface =5^ 
surrounding the object. 

In the following, we will still denote by Pexact the solution of ([2]), assuming 
that we know this latter for all x in and all t in R+. 
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2.3 Some useful facts about the standard wave equation 

In this section we give some basic properties of the following wave equation: 

Lou{x, t) = 0, (x, t) G X E+, 

u{x,0) = l{x), xeR^, (3) 

dtu{x,0) = h{x), xGR^, 

where I and h are two functions with compact support in B{0,1), 
and Lq := du — All this material can be found in (lit exam- 
ple. According to Huyghens' principle in the strong sense, for every cou- 
ple (xo,to) G IR-^ X the classical solution u{xQ,tQ) of ([3]) only depends 
on the values of / and h on the backward characteristic cone: 

||x — xoll = |t — tol- 

As a matter of fact, when the initial conditions / and h have their support 
in -6(0, 1), the solution at any time t has its support in -6(0, 1 + t). 

Huyghens' principle leads to the energy conservation for the solution n of ([3]). 

Definition 1. For every time T, we call energy of u the quantity: 

Eu{T) := ^ (lh(-,r) 11^1(5^3) + ||5tM(-,r)||^2(]a3)) . 

The variation of the energy of u is estimated from the following identity: 

3 

1=1 

Since u has a compact support, integrating for < t < T and over the whole 
space yields: 

Eu{T) = Eu{0). 
Note that if h is zero, we have: 

1 2 
Eu{0) = 2 IKIIh1(r3) • 



dtu {dttu - Au) = = dt 



\Vu\ 



+ 



{dtuf 
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2.4 Nudging 



Now, let us introduce the application of the BFN algorithm to the general 
TAT problem. 

Assume that the observation surface 5^ is included in the sphere S(0, 1). As 
exposed in Section 12. H our purpose is to compute an approximation of the 
original object /o from an incomplete set of data (/'Pdataj 

where pdata is a 

the solution (possibly contaminated by noise) of the wave equation ([2]) and 
where allows the knowledge of Pdata only on the observation surface 5^ . 
Theoretically speaking, should be 1^, but in practice the discretization 
process allows us to choose a C°° function with compact support. Indeed, 
given a resolution, there exists e > such that any function supported 
in r(=y,e) and equal to 1 on has the same discrete counterpart as 1^ 
(see Figure [1]). This is the reason why we will make the following assumptions 
on 0: 

i Vx e E^, < ^{x) < 1; 

ii (j) has a compact support in T(=5^,e); 

iii Vx e T{y,e/2), </)(x) = 1. 

Finally, in order to avoid some technical difficulties, we will assume that the 
initial object /o has his support in B(0, 1 — e) and is such that every Cauchy 
problem we shall consider has a classical solution, say C°°. 

We may proceed as follows: 

Definition 2 (Iterate of the BFN algorithm for TAT). 

1. [Forward evolution] 

From an rough estimate /o,i of the object to reconstruct Jq, defined as 
an approximate of the initial condition supported in B{0, 1— e), we first 
consider a solution pi of Problem ([3]) with initial conditions I = /o,i 
and h = (as usual in TAT). 

We compute this solution until time T = 2, chosen so that both pi{-,t) 
and Pexact(")i) (just as every function with compact support in the 
ball B(0, 1 — e)) vanish on B{0, 1 + e) for all t > T (cf Huyghens' 
principle in Section \2.3\) . 
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Figure 1: The body to investigate is surrounded by the acquisition surface 5^ ^ 
and generates a pressure wave. 



Thus at this final time, in the noise-free situation, the innovation 
term pi — Pdata? defined as the difference between the solution pi of 
Problem ([3]) and the data, is still a solution of ^ with I = fo — fo,i 
and h = and has an energy matching up to: 

1 2 

^Pl-Pdata(^) = -^Pl-Pdata(O) = 2 H^^-^'o ~ /o, 1 ) 1 1 1 (R3 ) : 

since the conservation of energy for the wave equation holds in this 
situation. 

2. [Backward evolution] 

Then we apply the backward nudging method to make the solution evolve 
back in time, with, as announced, the addition in the backward equation 
of a newtonian feedback, which adjusts the solution along its evolution 
by recalling it to the observed data. 

Namely, we add to L the feedback correction term k(j)dt{- — Pdata); fof 
some nudging parameter k G 11+ . And we compute the backward 
solution starting from the final state of the forward implementation, 
i.e. pi{x,T) and dtPi{x,T), untill the initial time t = 0. 

We obtain a corrected solution pi of the Cauchy problem, called back- 
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ward nudging problem: 

Lpi (X,t) = k(t){x)dt {pi - Pdata) (X,t), {x,t)eR^X [0, T] , 

pi{x,T) = pi{x,T), (4) 
dtPiix,T) = dtPi{x,T), xGR\ 

In order to obtain an equation in forward time, let us consider the 
map t ^ pi[T — t), still denoted by pi. The Cauchy problem becomes: 

Lpi = -k(t>{x)dt{pi+Pd.t.){x,T-t), {x,t) G R3 X [0,T], 
pi(x,0) =pi(x,r), x£R^ (5) 

dtpiix,0) = -dtPi{x,T), x£R^ 

where L is the operator L backwards in time. 

Obviously, in practice, when L is not the wave operator, there is no 
reason for this operator to be reversible. Nevertheless the newtonian 
feedback may act as a regularization term and keep the computation 
stable. 

3. [Update of the estimate] 

Finally, after this back and forth evolution, a new estimate pi{T) is 
obtained, but in order to iterate the process we need the new estimate 
to be supported in B{0, 1 — e), just like /o and /o,i. This is simply done 
by defining: 

fo,2{x) := tBio,i^e)Pi{x,T), X G E^. 

The last step can be easily justified by the fact that, since /o is supported 
in B{0,l-e): 

ll/o -pi(-,r)||^i(j^3) > ||/o - lB{0,l-e)Pl(-:^)||^l(]a3) • 

This scheme is then iterated, constructing a sequence of estimates (/o,n)nei:lN* ? 
which may converge to /q. The following theorem, proved in Section [3l en- 
sures that the convergence is geometrical in /7q(R^) under the assumptions: 

(i) L is the standard wave operator 9ft — A; 

(ii) The data are noise-free, i.e. Pdata = Pexact (see Section [2.ip . and the 
observation surface ^5^ is the whole sphere S{0, 1); 
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{in) The object to reconstruct /q is a C°° function with compact support 
ini?(0,l); 

(iv) /o,i and /0.2 are two consecutive estimates of /o obtained as described 
in Definition [2J 

Theorem 3. Under the assumptions {i)-{iv), there exists 1 > s > only 
depending on /q, (p the nudging parameter k > such that: 

II /O - /0,2||^^l(R3) <s\\fo- /0,i|Ih1{R3) • (6) 

We only consider this ideal framework for the sake of the proof, but numerical 
results given in Section [5] apply with more realistic assumptions and show a 
better convergence rate. 



3 Proof of Theorem [3] 

This section is devoted to the proof of Theorem [3j We study the case 
where L = Lq := — A is the standard wave operator. Let (p, /o and Pdata 
be as in Theorem [3j In particular, for the sake of the proof, we make the 
assumptions that the data are noise- free and that, since ,y = 5(0,1), the 
set T{y, e/2) is equal to ^4(1 - e/2, 1 + e/2). We shall use the notation: 

Lk:=du- ^ + k(f)dt, keR+. 

We first study the existence of the objects used in the BFN approximation 
scheme. 



3.1 Solution of the backward nudging equation 

The proof of the following theorem can be found in [7], Chapter 2.8. 

Theorem 4. Given three C°°-functions 0, / and h with compact support 
and k G 11+, T > there exists a unique C°° solution p of the Cauchy problem: 

Lkp{x,t) = -k(P(x)dtPd^t^{x,T -t), {x,t) G X [0,T], 

p{x,0) = l{x), x£R^, (7) 

dtp{x,0) = h{x), xeR^. 
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Moreover, for every (xo,to) G K.^ x p(xo,to) only depends on the values 
of /, h, (p and Pdata(")^ — ") in the backward characteristic cone: 

ll^; — 2;o|| < |t — tol- 

In particular, this means that for every sketch /o,i with compact support and 
of class C°°, the solution pi of the backward nudging equation ([7]) exists, and 
for all t £ [0,r] the function pi{-,t) has a compact support. This property 
will be helpful to apply an energy method. 

Now, it is easy to show that if pi is solution of ([7]) with initial condi- 
tions l{x) = pi{x,0) = pi{x,T) and h{x) = dtpi{x,0) = —dtPi{x,T), then 
the function u := pi — Pdatai',T — •) is solution of the Cauchy problem: 

Lku{x,t) = 0, {x,t) eR^ X [0,T], 

u{x,0) = pi{x,T) - Pda.ta.{x,T), X £ , (8) 

dtu{x,0) = -dt{pi - Pda.ta.){x,T), X £R^, 

and satisfies: 

1 2 
Eu{0) = Ep^^p^^^^{T) = - ||/o,i - /o||j:^i(R3) • 

Since, from the definition of /o,2: 

1 ~ 2 1 2 

Eu{T) > - ||pi(-,r) - /o||hi(r3) > 2 ll/o,2 - /o|Ihi(R3) ' 

we only need to prove that: 

Eu{T) < sEu{0), for some < s < 1, 

to complete the proof of Theorem [3j This is the purpose of the next section. 

3.2 Energy inequality 

We have the following identity for every solution n of ([7]): 

3 

+ k(t){dtuf - ^^dx.idtudx.u). 

i=l 

Recall that u has a compact support for every t, thus, integrating the latter 
equality over the whole space and for < t < T, we get: 

E^{T) - E^iO) = -k r f mnf- (9) 



dtuLkU = = dt 



\Vu\ 



+ 
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Consequently, the loss of energy of u is exactly the amount of its 'energy' 
passing through the support of during time T. 



In |12| |. in order to invert the spherical Radon transform, Finch et al proved 
the following trace identity, which will be usefull to estimate ([9]): 

Theorem 5. Suppose hi € C^{B{<d, p)) and Uj is the solution of ([3]) for Z = 
and h = hi, i = 1,2. Then we have the identity 



If,, -1 



oo 



„ , hih2 = — / / tuidttU2. 
^ Jr^ P Jo Js{q,p) 

From this we will obtain a key estimate of the 'energy' passing through the 
sphere 5(0,/)). Indeed, if u is solution of ([3]) with h = 0, taking hi = f 
and /i2 = A/ in the theorem yields: 

^ ^ fAf=— r [ t{dtu)\ 



2 Jrs p Jo Js(o,p) 

which gives: 



2 ' P Jo Js(0,p) 



'S(0,p) 

Let us introduce the solution uq of the Cauchy problem: 

Louo{x,t) = 0, {x,t) eR^ X [0,T], 

uo{x,0) = pi(x,T) -pdata(a;,r), x G R^, 

uot{x,0) = -dtPi{x,T) + (?iPdata(a;,0), X e R^, 

which is nothing but {pi - Pdata)(-,7' - •). 

Then u = uq + v where v is solution of: 

Lov{x,t) = -k(j)dt{uo + v){x,t), {x,t) £ R-^ x [0,T], 
v{x, 0) = dtv{x, 0) = 0, X G R^. 

The calculation of the energy of v yields: 

E,{T) = -k[ [ ^dtuodtv-k [ [ cp{dtvf. 
Jo JR3 Jo Jr3 

Thanks to the Cauchy-Schwarz inequality, one has: 

E.{T)<k y^dtuo ^^^^^^^ H{T)-kH{Tf, (11) 
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where: 



Moreover, we have: 



H{T):= ^dt' 



L2([o,T]x]R3) 



E.{T)>f {dtvf{;T)> [ ct>{dtvf{;T) = ^{H^){T). 

J]R3 J^A dT 

So that, dividing by H{T) in ([TT]) we find: 

4^HiT) + kHiT) < k J^dtuo 
A classical calculation shows that this diff'erential inequation yields: 











Jo 





L2([0,s]xR3) 



and since the map s i— )• || V^f^t^o ||^2([q ^jxrs) increasing with s, we finally 
get: 

This inequality expresses the fact that despite the feedback part of Equa- 
tion ([7]), its solution u keeps a certain amount of energy on the support of (j), 
which is proportional to the same energy for uq. Indeed: 



L2([0,r]xR3) 



> 



> 



L2{[0,T]xR3) 



L2([0,r]xR3) 



H{T) 



-kT 



thus: 



JR3 



JR3 



(t){dtUo) e 



2 \ -2kT 



(12) 



The last step of this proof is the estimation of A. We have: 
A> [ r{dtUof= [ r{dtPi-dtPd.t.f{T-t), 

J A(l-e/2,l+e/2) JO J A(l-e/2,l+e/2) JO 



IA{l-£/2,l+s/2) 
SO that: 



A > 



IA(l-£/2,l+e/2) JO 



{dtpi-dtPda.ti>.y 



l+e/2 pT 



l-e/2 JO JS{0,p) 



{dtPl- dtp Aa.taf ■ 
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Since pi — Pdata vanishes on ^(1 — e/2, 1 + e/2) for every t > T: 

rl+e/2 r-co r- 

A>- / / t{dtpi - dtPd^t.?. 



l-e/2 Jo JS{0,p) 



T 

Finally, applying ^Qii to pi - pdata, one has: 

1 I 



1 1 2 

which yields: 



A > ^Eu{0). (13) 
Combining the above estimate with ([9]) and (|12p . we get: 

Eu{T) < (1 - ^e'''^)E^{0), 

which was the inequality to prove, with s := 1 — ^e~'^^'^ < 1. 

Remark 6. Since we know that T = 2, there exists an optimal choice for the 
nudging parameter k which leads to s = 1 — Even though this theoretical 
convergence rate does not seem satisfactory, we shall see in Section that 
the actual convergence is much faster in practice. 



4 Variational approach 

We deal now with a variational formulation of the inverse TAT problem, that 
leads to some new reconstruction techniques. 

In this section, we consider the reconstruction of an object /q G Hq{B(0, 1)) 
from a set of observations pdata = Pexact +Pnoise- Here, pexact := ^/o where 
the linear operator maps an initial condition / to its related solution of the 
Cauchy problem ^ with a null derivative initial condition h = 0, restricted 
to a non-empty, open observation set C A{1 — e,l + e). Namely, W can 
be written W = ^W, with: 

W: Fi(i?(0,l)) C'{[0,T];H^{B{0,R))) 
/ ^ Wf, 

where Wf is the weak solution of ([3]) with the initial conditions I = f 
and h = 0. Even though Wf is defined on the whole space R^, we can 
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choose R > large enough for any solution of ([3]) to keep a null trace 
on S{0,R), such that the restriction of Wf to B{0,R) is in H^{B{0,R)) for 
every t £ [0, T]. The operator $ is the restriction operator: 

^: L^{[0,T]xR^) l2([0,T]x^,) 
9 ^ 5U,- 

Our purpose is to solve the inverse problem: 

= Pdata, 

by means of the following minimization problem: 

Minimize Ja{f) 

s.t. f eHl{B{Q,l)). 

Here, is the Tikhonov functional: 

J«: ^1 (5(0,1)) E+ 

1 f'^ 2 OL 2 

f ' — ^2/ "'^•^~^^'^*^"^''(^^) 2" ' 



with a > 0. Since is a bounded linear operator, is a strictly con- 
vex (thanks to the regularization parameter a), differentiable and coercive 
functional, and any solution of {^a) is characterized by: 

VJ«(/) = W*{Wf-pd,t.) + af = 0. 
For any a > 0, the above equation has a unique solution: 

so that Problem {^a) is well-posed. We have chosen to solve this latter by 
means of the conjugate gradient method, which requires, at each iteration, 
the computation of VJ^. 

Now, one has, for every i; G C^{B{0, 1)) and / e H°{B{0, 1)): 

(VJ,(/),^)i2(B(0,l)) = (^*(#^/-pdata),^)+a(/,V) 
= (^/-Pdata,^V')+a(/,V'), 
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with: 



G 



T 



J 

T 



JB{0,R) 



lA,(t^/-Pdata)VFV'. 



Considering the weak solution u* of the adjoint equation: 



dttu*ix,t) - Au*{x,t) = lAAWf-Pd^t^), 

n*(x,T) = 0, (14) 

dtu*{x,T) = 0, 



and since Wip G Cq°{B{0,R)), we have: 

Gf,^= [ [ dttu*Wi^+ [ [ (Vn*,VVFV). 
Jo Jb{o,r) Jo Jb{o,r) 

Since B{0,R) has been chosen large enough, applying Green's formula and 
integrating by parts yield: 



G 



Finally: 



u*LWij- / dtu*{x,0)Wi^{0)dx 
Jb{o,r) Jb{o,r) 

{-dtU*{;0),,P). 



VJM) = -dtu*i;0) + af. 



This calculation ensures that each iteration of the conjugate gradient algo- 
rithm necessitates successively the computations of Wf and of its adjoint 
state u* , whose algorithmic complexity is comparable to one iteration of the 
BFN algorithm. This should be noted that the computation of the adjoint 
state requires the storage of Wf — Pdata on A^, which is not required in the 
BFN algorithm. 

Note that even if we have formulated the calculation of V Jq in a continuous 
framework, its discrete counterpart leads to similar considerations. Obvi- 
ously, the behaviour of the method strongly depends on the scheme chosen 
to compute solutions of i^. 
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5 Numerical implementation 



In order to test the numerical behavior of the BFN method, we have imple- 
mented it in several situations. In particular, we have considered the case 
where the data are available only on a half-sphere, so that the illustration 
fits the practical experimental set-up. 



5.1 Implementation framework 

We now describe the numerical context chosen to compare Back and Forth 
Nudging (BFN) and the Conjugate Gradient (CG) techniques to a reasonable 



algorithm, namely the Neumann Series (NS) (since it is shown in [33| that 
the use of NS significantly improves the Time Reversal (TR) algorithm, we 
mainly chose to consider NS results). 

We work now in a 2D framework for numerical costs reasons. The object to 
reconstruct is in [— 0.5,0.5]'^ gridded with 256^ pixels, such that the space 
step is 5x = 1/256. 

Sensors are located on the circle S{0, surrounding the object, on which 

various distributions are considered. If no additional information is given, 
800 detectors are used when the data are complete. Nevertheless some re- 
constructions were computed with less detectors in order to bring to light the 
robustness of the techniques. When a multiplicative white Gaussian noise is 
added to the data, their level is set to 15%. 

The final time T = y/2 is then taken minimal in order to leave waves traveling 
to the sensors when sound speed is constant with value 1. Note that the four 
studied methods can significantly be improved by increasing T (c/. [33]). We 
use the domain [—0.5 — \/2 — e; 0.5 + -v/2 + e]^ for the computation, where e 
represents a few pixels, which is large enough to avoid refiection effects on 
its boundary. 

The model is computed by means of the classical finite difference time domain 
(FDTD) method (explicit Euler scheme and classical five points discretiza- 
tion of the Laplace operator A^). 

Figure [2] represents the 2-dimensional objects to reconstruct, latter denoted 
by /q. Squares able to give rise to typical phenomenons of observation of 
wave-like systems while the second object, known as the Shepp- Logan phan- 
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torn, allows good analysis and comparisons of the ability of the methods. 
The third object has been considered for combining sharp contom's with 
smoother areas. Each is implemented as a 256 by 256 matrix. 




Figm'e 2: The initial objects /q: squares, Shepp-Logan and skull. 



Four sound speeds are used. First a constant one c = 1, then three other 
one, variable, which are chosen as in [s^, respectively defined as: 

c{x,y) = 1 + 0.2sin(27ra;) + 0.1cos(27ry), 
1 + 9(x"^ + y^j 

-0.4exp(-10(3\/x2 + y2 _ 2)2), 
c{x,y) = 1.25 + sin(27r2;) cos(27ry). 

They are plotted on figure [3] and abbreviated respectively NTS, TSl and 
TS2, as the first is a Non- Trapping Speed while the two other are Trapping 
Speeds. No cutoff is used to smooth the transition between variable internal 
speed and constant external speed (always being 1). 




Figure 3: The different variable speeds: NTS, TSl and TS2 (from equa- 
tions (fTSjl 'l. 
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Concerning the parameters of the algorithms, the regularization parameter a 
in CG method equals 0. Other values were tested that stabilized the algo- 
rithm but always blur the reconstructions so that we get higher errors than 
without regularization. The observed convergence rates of CG and NS meth- 
ods are fairly high, so that only 10 iterations yield a good approximation of 
the best achievable solutions. With regard to the two nudging coefficients, 
while we can theoretically affirm that the larger it is, the faster the solution 
tends to the observations (thanks to the energy equality the discretiza- 
tion of the time-derivative of the solution imposes a coefficient k not larger 
than l/6t. In order to obtain the best convergence and to make the choice 
of k less arbitrary (and fitting with the limits of the numerical method), it 
is defined as 0.9/ 6t. 

The following tables give the relative mean square error obtained with the 
best estimation obtained after at most 10 iterations and, in brackets, the 
corresponding number of iterations required. Note that BFN, CG and NS 
furnishes convex error plots with a minimum. 



Finally, as it has been suggested in [3j|, we also have used an artificial 
numerical attenuation to offset spurious high-frequency effects and noise, 
with interesting consequences. Then for BFN, NS and CG techniques, the 
numerical model becomes: 

^ = AsxUn + OX^Afe — , 

where the value for e is empirically set and an optimal value appears for each 
method, which are approximately 2 for BFN and 1.7 for CG and NS in a 
complete data situation. 

The attenuation being only introduced in a numerical improvement view, it 
is set to attenuate both back and forth solutions, should the case arise, and 
when the data are corrupted by noise. 



5.2 Implementation results 

The first results are given in Table [T] in a noiseless complete data situation. 
One sees that variable speeds can have various effects, depending on the 
method we use. NS offers the best reconstruction, then BFN and finally 
CG. 

These techniques can be consequently improved by increasing the number of 
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Speed map BFN 

c=l 2.3(10) 

NTS 3.2(10) 

TSl 7.1(10) 
TS2 7(10) 



CG 
7.1(3) 
7.5(4) 
16.5(5) 
10.7(5) 



NS 
1.4(10) 
2.3(10) 
1.6(10) 
6.6(10) 



Table 1: Noiseless complete data for Shepp-Logan phantom: relative mean 
square error (number of iterations, limited to 10). 

iterations as we can see in Table [2] for the BFN, where 100 iterations instead 
of 10 are computed. 



Table 2: Noiseless complete data from Shepp-Logan phantom: relative mean 
square error (BFN, 100 iterations). 

When adding a 15% level noise in the constant speed situation, one sees on 
Table [3] that any method offers an interesting reconstruction: one gets less 
than 20% error for 15% noise thanks to the attenuating scheme. 



Table 3: Noisy complete data (15% noise) from Shepp-Logan phantom with 
the constant speed. Reconstruction with attenuating sheme: relative mean 
square error (number of iterations, limited to 10). 

About the propagation of the singularities and their observation, we put 
the light on the effect of getting incomplete data (in a noiseless situation) 
in Fig. m The sensors are located on the upper half circle, that illustrate 
the practical case where the observed body cannot be surrounded by the 
observation surface. 

In this situation (constant speed and noiseless data on the half circle), the 
three methods are compared in Fig. [5] and [6l NS is excellent, while BFN 



c=l NTS TSl TS2 
1 1.5(100) 1.2(100) 4.1(100) 



BFN CG NS 
17.6(10) 19.8(2) 17.1(2) 



21 



Figure 4: Reconstructions by 100 BFN iterations. Squares observed singu- 
larities. Complete (1.1% error) and half-circle (10.4%) noiseless data with 
constant speed. 

equivalent to CG convergence in about 5 iterations, but keeps on improving 
the reconstruction along next iterations. As always, even if its reconstruction 
is the worst one, CG converges in less iterations and yields a good rough 
estimate faster than the others. 




Figure 5: Incomplete noiseless data on the half circle from Shepp-Logan 
phantom with constant speed: relative mean square error plots (BFN in 
solid line, CG in dashed-dotted line and NS in dashed line). 

Still testing the robustness of the methods, we reduce the number of sensors. 
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BFN Reccnjtiudion 



TR Reconitruciion 




BFN CG TR NS 
10.6(10) 16(6) 52 6.2(10) 



Figure 6: Incomplete noiseless data on the half circle from Shepp-Logan 
phantom with constant speed: reconstructions with BFN (top left), TR (top 
right), CG (bottom left) and NS (bottom right) and relative mean square 
errors table (error, number of iterations, limited to 10). 

which are still homogeneously located on the observation circle, but only one 
out of a is considered in the data. The results shown in Table 3] highlight 
some redundancy aspect of the data since the reconstructions have a similar 
quality in all situations, particularly with CG which is very robust here. 
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That is obvious at convergence of the algorithms, but one must consider 
more than 100 iterations for BFN and NS when a = 8. Finally, if the 
intensity is not well recovered, forms are reconstructed with few artifacts. 



Settings 


BFN 


CG 


NS 


a = 1 


2.3(10) 


7.1(3) 


1.4(10) 


a = 2 


5.4(10) 


8.9(4) 


4.5(10) 


(7 = 4 


12.1(10) 


12.4(4) 


11.5(10) 


a = 8 


25.3(10) 


15(5) 


18.5(10) 



Table 4: Partial noiseless data (one out of a sensor records data) from Shepp- 
Logan phantom with constant speed: relative mean square errors (number 
of iterations, limited to 10). 

Then we sum up the different unfavorable conditions from above and show 
the results on Fig. [7] and [8j Here data are acquired on the half circle, only 
one out of two sensors is considered, a 15% level noise is added and the speed 
is constant. The last object, a skull, has been chosen in this last experiment 
in order to consider a more complex structured object. 




Figure 7: Incomplete partial data {a = 2) on the half circle with 15% noise 
from skull phantom and with constant speed: relative mean square error 
plots (BFN in solid line, CG in dashed-dotted line and NS in dashed line). 

Finally, even if the reconstructions obtained in this last situation are not 
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BFH Reconstiodion 



TR Reconjiiuciion 




BFN CG TR NS 
30.7(4) 37.8(10) 59.9 37.9(4) 

Figure 8: Incomplete 15% noise data on the half circle from skull phantom 
with constant speed and partial data [a = 2): reconstructions with BFN 
(top left), TR (top right), CG (bottom left) and NS (bottom right) and 
relative mean square errors table (error, number of iterations, limited to 10). 

satisfying, it may fit with very bad real conditions (where the noise level can 
be far greater than 15% for example) and allows us to test the robustness of 
the methods. In this view. Table [5] shows the errors if we consider variable 
speeds in a similar case but with a = 1. 
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Speed map BFN 

c = 1 36.4(8) 

NTS 36.3(8) 

TSl 46.3(4) 

TS2 39.1(8) 



CG 
40(10) 
40.6(10) 
50.9(8) 
42.3(9) 



NS 
38.7(3) 
38.6(3) 
51.9(3) 
41.5(4) 



Table 5: Incomplete 15% noise data on the half circle from Shepp-Logan 
phantom: relative mean square errors (number of iterations, limited to 10). 

5.3 Implementation conclusions 

CG reacts much badly to noise addition above all, but it provides good re- 
constructions after a couple of iterations in comparison with BFN and NS, 
so that it can be used to get a rough estimate for another method. More- 
over CG keeps interesting when data are incomplete, not overreacting. On 
the whole, BFN and NS provide good and quite equivalent reconstructions 
with a similar cost. Finally, NS furnishes better reconstruction, but appears 
less robust than BFN to noise and to restriction of the number of sensors. 
Numerical attenuation provides a good regularization that offsets noise for 
NS and BFN. Notice that this artificial numerical attenuation may allow us 
to consider lossy medium in back and forth evolutions as the physical loss 
does not exceed the artificial one. 

Depending on whether data are complete or not, one may choose CG to get 
a first rough estimate, then NS or BFN to reconstruct the object. It has to 
be said that not stopping iterations to 10 often improves results in noiseless 
and/or complete data situations, but that this number of iterations provides 
the best estimate in less favorable settings. 

6 Conclusions 

We have established the theoretical geometrical convergence of the TAT- 
BFN algorithm, with an explicit convergence rate in an ideal situation. In 
principle, this method enables to make weak assumptions in the TAT model- 
ing (variable sound speed, damping, incomplete data, etc.), since the model 
is considered as a weak constraint. Even though the proof of the convergence 
has only been stated in an ideal framework, the numerical implementation 
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of the algorithm showed that the BFN method remains efficient in more 
general situations. Indeed this technique still yields robust reconstructions 
in the noisy and incomplete data situation. Even though the convergence 
rates of the Conjugate Gradient and of the Neumann Series methods ap- 
pear to be better, the BFN provides more satisfactory asymptotic behavior 
and minimal relative error, along with a low numerical complexity, in quite 
unfavorable observation situations. 

From these preliminary results, one should now focus on their generaliza- 
tion to more realistic wave equations. For example, one may introduce some 
damping in the model and study the theoretical validity and convergence of 
the methods. Our first attempts in this direction (mainly numerical imple- 
mentations) suggest future interesting developments. 
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